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1  SUMMARY 


Contributions  of  this  research  are  the  development  of  an  on-board  relative  motion  maneuver  plan¬ 
ning  approach  for  spacecraft  debris  avoidance  that  can  handle  set  bounded  disturbances.  In  this 
report  we  describe  the  development  of  an  on-board  maneuver  planning  approach  based  on  the 
use  of  constraint-admissible  positively  invariant  sets  [1].  The  sets  determine  connectivity  between 
forced  and  unforced  spacecraft  equilibria  in  the  Clohessy-Wiltshire-Hill  (CWH)  relative  motion 
frame  [2],  The  collection  of  equilibria  form  a  virtual  net  in  the  vicinity  of  the  spacecraft.  Two 
equilibria  are  connected  if  a  choice  of  a  Linear  Quadratic  (LQ)  feedback  gain  can  be  made  that 
results  in  a  transition  between  the  equilibria  which  avoids  collision  with  a  potentially  moving  de¬ 
bris/obstacle  while  satisfying  limits  on  thrust.  A  connectivity  graph  for  all  the  equilibria  in  the  net 
is  constructed  based  on  fast  growth  distance  computation  between  two  ellipsoidal  sets,  while  real¬ 
time  graph  search  algorithms  are  used  to  optimize  an  equilibria  hopping  sequence  to  avoid  debris 
collisions.  Unlike  existing  spacecraft  trajectory  optimization  techniques,  our  method  does  not  rely 
on  precise  assignment  of  spacecraft  position  along  the  trajectory,  and  is  able  to  assure  robustness 
to  unmeasured  (but  set-bounded)  disturbances  and  uncertainties. 

The  following  papers  have  been  published  related  to  the  subject  matter  of  this  report. 

•  A.  Weiss,  M.  Baldwin,  R.  S.  Erwin,  and  I.  Kolmanovsky,  “Spacecraft  Constrained  Maneuver 
Planning  for  Moving  Obstacle  Avoidance  Using  Positively  Invariant  Constraint  Admissible 
Sets,”  Proc.  Amer.  Contr.  Conf.,  Washington,  DC,  June  2013. 

•  A.  Weiss,  M.  Baldwin,  R.  S.  Erwin,  and  I.  Kolmanovsky,  “Spacecraft  Constrained  Maneu¬ 
ver  Planning  Using  Positively  Invariant  Constraint  Admissible  Sets,”  Bar-Itzhack  Memorial 
Symposium,  Haifa,  Israel,  October  2012. 

•  M.  Baldwin,  A.  Weiss,  I.  Kolmanovsky,  and  R.  S.  Erwin,  “Spacecraft  Debris  Avoidance 
using  Positively  Invariant  Constraint  Admissible  Sets,”  AAS/AIAA  Space  Flight  Mechanics 
Meeting,  Charleston,  SC,  January  2012,  AAS- 12-250  . 

The  project  has  supported  one  doctoral  student  (Avishai  Weiss)  for  one  year  of  his  doctoral 
studies. 


2  INTRODUCTION 

Orbital  debris  is  a  growing  problem,  with  about  40%  of  ground-trackable  objects  originating  from 
explosions  that  now  number  approximately  5  per  year  [3].  Spacecraft  maneuver  planning  proce¬ 
dures  thus  have  to  address  debris  avoidance  requirements.  While  obstacle  avoidance  is  a  standard 
problem  in  robotics  [4,5],  the  related  spacecraft  problems  have  several  unique  features.  In  partic¬ 
ular,  the  space  environment  is  relatively  uncluttered,  thus  permitting  for  a  variety  of  maneuvers. 
Spacecraft  dynamics  are  quite  different  from  those  of  typical  robots.  Maneuver  efficiency  with 
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respect  to  time  and  fuel  consumption  is  a  critical  consideration.  The  states  of  the  spacecraft  and 
the  debris  can  only  be  estimated,  often  with  a  significant  estimation  error.  Finally,  computational 
algorithms  must  be  fast  and  optimized  given  moving  objects  and  the  limited  computing  power  on¬ 
board  most  spacecraft.  These  unique  features  of  spacecraft  maneuver  planning  problems  provide 
the  motivation  for  the  development  of  specialized  algorithms. 

Interest  in  spacecraft  trajectory  optimization  with  obstacle  avoidance  has  increased  in  recent 
years.  An  optimal  control  problem  with  path  constraints  constructed  as  ‘keep  out’  zones  to  avoid 
obstacles  was  formulated  in  [6].  The  Sparse  Optimal  Control  Software  (SOCS)  software  was  then 
used  to  solve  the  problem  [7].  Another  nonlinear  optimal  control  formulation  was  used  in  [8] 
to  solve  for  minimum-fuel  rendezvous  between  a  target  and  chaser,  where  collision  avoidance 
requirements  were  incorporated  as  inequality  constraints.  The  method  involved  solving  a  sequence 
of  unconstrained  optimal  control  problems,  whose  solution  converges  to  the  solution  of  the  original 
problem.  A  3-D  static  optimization  over  final  relative  position  and  time-of-flight  such  that  obstacles 
are  avoided  and  cost  is  optimized  is  presented  in  [9].  Feedback  is  incorporated  by  re-planning  over 
either  constant  or  variable  time  intervals. 

Debris  avoidance  strategies  have  also  been  defined  utilizing  collision  avoidance  probabilities. 
Collision  avoidance  strategies  based  upon  the  number  of  evasive  maneuvers,  expected  risk  reduc¬ 
tion,  false  alarm  rate,  required  propellant  consumption,  and  mass  fraction  for  an  accepted  collision 
probability  are  presented  in  [10]. 

Guidance  based  on  artificial  potential  function  is  used  in  [9, 1 1]  to  determine  a  rendezvous  path 
free  of  obstacles.  A  potential  function  is  developed  with  the  intent  that  a  minimum  occurs  at  a 
desired  relative  position  and  then  a  dynamic  control  law  is  used  to  ensure  the  trajectory  is  obstacle 
free  [11]. 

The  spacecraft  obstacle  avoidance  problem  has  also  been  treated  using  linear  programming 
techniques  [12-15].  In  [12],  the  minimum-fuel  avoidance  maneuver  is  formulated  with  linear 
constraints  and  discrete  dynamics  modeled  as  an  linear  time-varying  (LTV)  system.  In  [13],  the 
trajectory  optimization  problem  is  formulated  as  a  linear  programming  problem  with  the  capability 
of  including  operational  constraints  and  the  optimal  number  of  maneuvers  is  determined.  In  [14], 
a  mixed-integer  linear  program  results  from  combining  collision  avoidance,  trajectory  optimiza¬ 
tion,  and  fleet  assignment  to  obtain  the  optimal  solution  for  spacecraft  maneuvers.  A  robust  linear 
programming  technique  is  proposed  in  [15].  The  maneuver  can  be  constructed  by  solving  a  lin¬ 
ear  programming  problem  with  no  integer  constraints  and  guaranteeing  collision  avoidance  with 
respect  to  bounded  navigation  uncertainty. 


3  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

We  develop  a  chained  invariant  set  method  to  avoid  both  static  and  moving  debris  during  spacecraft 
relative  motion  maneuvers.  The  equations  of  motion  for  spacecraft  orbital  dynamics  are  reviewed 
in  Appendix  A.  We  use  linearized  equations  of  motion  (A.3),  where  continuous  thrust  actuation 
Uk  =  [Fx,k,  FyJc,  FzM]T  is  assumed. 
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Our  approach  to  debris  avoidance  is  based  on  utilizing  constraint-admissible  positively  invari¬ 
ant  sets  [1, 16-18]  centered  around  spacecraft  forced  and  unforced  equilibria.  A  finite  set  of  these 
equilibria  used  for  constructing  debris  avoidance  maneuvers  is  referred  to  as  a  virtual  net.  Given 
an  estimate  of  the  debris  position,  we  build  a  connectivity  graph  that  identifies  the  equilibria  in 
the  virtual  net  between  which  the  spacecraft  can  move,  with  guaranteed  collision-free  motion  and 
within  the  available  thrust  authority.  We  then  employ  graph  search  to  determine  an  efficient  path 
between  the  equilibria  that  ensures  debris  avoidance.  One  of  the  main  reasons  this  framework  is  at¬ 
tractive  compared  to  alternatives  such  as  open-loop  trajectory  planning,  is  the  ability  to  incorporate 
bounded  disturbances  such  as  thrust  errors,  air  drag,  and  solar  pressure. 

3.1  Virtual  Net 

The  virtual  net  comprises  a  finite  set  of  equilibria,  Xe(r),  corresponding  to  a  finite  set  of  prescribed 
spacecraft  positions  r  e  J\f  =  {ri ,  ri, . . . ,  rn)  c  R3 , 

Xe(rk)  =  \rk  0  ]T  =  [  fx,k  ry,k  rz,k  0  0  0  ]T,  k  =  1, •  •  •  ,n,  (1) 

whose  velocity  states  are  zero,  and  where  n  is  the  number  of  equilibria  in  the  virtual  net.  See 
Figure  1.  We  assume  that  for  all  rejV,  the  corresponding  values  of  control  necessary  to  support 
the  specified  equilibria  in  steady-state  satisfy  the  imposed  thrust  limits. 


Figure  1:  The  virtual  net  for  debris  avoidance.  Dots  correspond  to  positions  at  equilibria,  Xe(r),  on 
a  virtual  net.  The  ellipsoid  represents  the  debris  position  and  uncertainty. 


3.2  LQ  Controller  with  Gain  Switching 

A  conventional  Linear-Quadratic  (LQ)  feedback 

U  =  K(X  -  Xe(r))  +  Yr  =  KX+  H(K)r,  (2) 
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is  used  to  control  the  spacecraft  thrust  to  arrive  at  a  commanded  equilibrium  (1),  where 


r  = 


-3 n2mc  0  0 

0  0  0 
0  0  n2mc 


h(K)  =  t-k  ^  , 

and  where  I3  denotes  the  3x3  identity  matrix  while  O3  denotes  the  3x3  zero  matrix.  The  LQ  con¬ 
troller  provides  an  asymptotically  stable  closed-loop  system  but  does  not  enforce  debris  avoidance 
constraints. 

To  provide  greater  flexibility  in  handling  constraints,  a  multimode  controller  architecture  is 
employed  [16].  Specifically,  we  assume  that  a  finite  set  of  LQ  gains  K  e  K,  =  , Km }  is 

available  to  control  the  spacecraft.  By  using  a  large  control  weight  in  the  LQ  cost  functional, 
motions  with  low  fuel  consumption  yet  large  excursions  can  be  generated;  using  a  large  control 
weight  in  the  LQ  cost,  motions  with  short  transition  time  can  be  generated  [19].  We  assume 
that  a  preference  ordering  has  been  defined  and  the  gains  are  arranged  in  the  order  of  descending 
preference,  from  K\  being  the  highest  preference  gain  to  Km  being  the  lowest  preference  gain. 


3.3  Positively  Invariant  Sets 

The  ellipsoidal  set 


C(r,K )  =  [X  e  R6  :  l-(X  - Xe(r))T P(K)(X -  Xe(r))  <  1}  c  R6, 


(3) 


where 


A(K)T PA(K)  -  P  <  0,  (4) 

A(K)  =  (A  +  BK ),  and  P  =  P(K)  >  0  is  positively  invariant  for  the  closed-loop  dynamics.  Positive 
invariance  implies  that  any  trajectory  of  the  closed-loop  system  that  starts  in  C(r,  K)  is  guaranteed 
to  stay  in  C(r,  K)  as  long  as  the  same  LQ  gain  K  is  used  and  the  set-point  command  r  is  maintained. 
To  achieve  the  positive  invariance,  the  matrix  P  can  be  obtained  as  the  solution  of  the  discrete-time 
Riccati  equation  in  the  LQ  problem  or  as  the  solution  of  the  above  Lyapunov  equation  for  the 
closed-loop  asymptotically  stable  system.  We  note  that,  because  the  system  is  linear,  the  positive 
invariance  of  C(r,K )  implies  the  positive  invariance  of  the  scaled  set 

C(r,  K,p )  =  {XeR6:l-{X-  Xe{r))T  P(K)(X  -  Xe(r))  <  p\  p  >  0. 
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Geometrically,  the  set  C(r,  K,p)  corresponds  to  an  ellipsoid  scaled  by  the  value  of  p  and  centered 
around  Xe(r),  r  ej\f. 


3.4  Debris  Representation 

We  use  a  set,  0(z,  0,  centered  around  the  position  z  e  R3,  to  over-bound  the  position  of  the  debris, 
i.e., 


0(z,Q)  =  {X  e  R6  :  (SX-zfQ(SX-z)  <  1}, 


(5) 


where  Q  =  QT  >  0  and 


S  = 


1  0  0  0  0  0 
0  1  0  0  0  0 
0  0  1  0  0  0 


(6) 


The  set  0(z,  Q )  can  account  for  the  debris  and  spacecraft  physical  sizes  and  also  for  the  uncertain¬ 
ties  in  the  estimation  of  the  debris/spacecraft  position.  Note  that  the  set  0(z,  Q )  has  an  ellipsoidal 
shape  in  the  position  directions  and  it  is  unbounded  in  the  velocity  directions.  Ellipsoidal  sets, 
rather  than  polyhedral  sets,  are  used  here  to  over-bound  the  debris,  since  ellipsoidal  bounds  are 
typically  produced  by  position  estimation  algorithms,  such  as  the  Extended  Kalman  Filter  (EKF). 


3.5  Static  Debris  Avoidance  Approach 

Consider  now  r,-  e  Af,  representing  a  possible  position  on  the  net  that  the  spacecraft  can  move  to  as 
a  part  of  the  debris  avoidance  maneuver.  Suppose  that  the  current  state  of  the  spacecraft  is  Z(to)  at 
the  time  instant  to  e  Z+ .  If  there  exists  a  p  >  0  and  Kj  e  K,  such  that 


X(t0)  e  C(n,  Kj,p )  and  0(z,  0  n  C(n,  Khp)  =  0,  (7) 

the  spacecraft  can  move  to  the  position  r,  €  M  by  engaging  the  control  law  with  r(t)  =  r,  and  K(t)  = 
Kj,  t  >  to,  and  without  hitting  the  debris  confined  to  0(z,  0.  This  idea  underlies  our  subsequent 
approach  to  debris  avoidance,  where  we  maintain  the  spacecraft  within  a  tube  formed  by  positively 
invariant  sets  that  do  not  intersect  with  debris. 

The  minimum  value  of  p  >  0  for  which  0(z, Q)f]C(r,K,p)  /  0  is  referred  to  as  the  growth 
distance  [20].  This  growth  distance  can  also  be  viewed  as  the  least  upper  bound  on  the  values  of 
p  for  which  0(z,  0  and  C(r,K,p )  do  not  intersect.  See  Figure  2.  We  use  the  notation  pg(r,K,  Q,z ) 
to  reflect  the  dependence  of  the  growth  distance  on  the  set-point  reAf,  the  control  gain  K  eJC  and 
the  obstacle  parameters  Q  and  z. 

Note  that  the  growth  distance  depends  on  the  position  of  the  debris  which  may  be  unknown  in 
advance.  Consequently,  growth  distance  computations  have  to  be  performed  online. 

Since  spacecraft  have  limited  thrust,  we  additionally  define  a  maximum  value  of  p  =  pu(r,K) 
for  which  X  e  C(r,  K,pu(r,  K))  implies  that  the  thrust  U  =  KX  +  H(K)r  satisfies  the  imposed  thrust 
limits.  We  refer  to  pu  as  the  thrust  limit  on  growth  distance.  Unlike  pg,  the  value  of  pu  does  not 
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Figure  2:  The  positively  invariant  set  is  grown  till  it  touches  the  debris. 

depend  on  the  position  or  shape  of  the  debris  and  can  be  pre-computed  off-line. 

Finally,  we  define  the  thrust  limited  growth  distance 

P*(r,  K,  Q,z )  =  min  {pg(r,  K,  Q,z),pu(r,  K)}.  (8) 

Note  that  X(to )  e  C(ri,Kj,p*(ri,Kj,z ))  implies  that  the  ensuing  closed-loop  spacecraft  trajectory 
under  the  control  (2),  where  r(t)  =  ri  and  K(t)  =  Kj  for  t  >  to,  satisfies  the  thrust  limits  and  avoids 
collisions  with  a  debris  confined  to  0(z,  0. 

The  above  definitions  were  given  for  the  case  of  a  single  stationary  debris,  0(z,  0.  In  the 
case  of  multiple  debris,  the  growth  distance  is  replaced  by  the  multi-growth  distance,  which  is  the 
minimum  growth  distance  to  each  of  0(zu  Qi ),  /  =  1, •  •  •  ,n^. 

3.6  Growth  Distance  Computations 

Define  X  =  X  -  Xe(r)  and  a  =  2 p2.  The  problem  of  determining  the  growth  distance  pg(r,  K,  Q,  z), 
reduces  to  the  following  constrained  optimization  problem: 


min  a 

a,X 

subject  to  XTPX  <  a  (9) 

((5  (X  +  Xe(r))  -  z)TQ((S (X  +  Xe{r))  -  z)  <  1 

To  solve  this  optimization  problem,  we  use  the  Karush-Kuhn-Tucker  (KKT)  conditions  [21,22], 
Note  that  standard  linear  independence  constraint  qualification  conditions  hold  given  that  P  >  0. 
We  define 


C  =  a  +  AxOPPX -  a)  +  A2((S (X  +  Xe(r))  - z)TQ(S (X  +  Xe(r))  -  z)  - 1), 
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where  A\  and  A2  are  Lagrange  multipliers.  The  stationarity  of  the  Lagrangian  (setting  partial 
derivative  equal  to  zero)  with  respect  to  a  yields  A  \  =  1 .  The  stationarity  of  the  Lagrangian  with 
respect  to  X  yields 


X  =  X(A2,r,z)  =  -(P  +  A2STQSrlSTQ(SXe(r)-z)A2,  (10) 

where  the  scalar  A2  >  0  is  to  be  determined.  Note  that  P  >  0,  STQS  >  0,  A2  >  0  (as  the  Lagrange 
multiplier  corresponding  to  an  inequality  constraint)  imply  that  (P  +  A2STQS)  is  invertible.  The 
problem  reduces  to  finding  a  nonnegative  scalar  A2,  which  is  the  root  of 

F(A2,r,z)  =  ((SX-z)tQ(SX-z)-  1  =  0,  (11) 


where 

X  =  X(A2,r,z)  +  Xe(r). 

The  scalar  root  finding  problem  (11)  has  to  be  solved  online  multiple  times  for  different  rej\f, 
and  in  the  case  of  avoiding  a  predicted  debris  path  also  for  different  z’s.  To  solve  this  problem 
fast,  while  reusing  previously  found  solutions  as  approximations,  a  dynamic  Newton-Raphson’s 
algorithm  is  used  [22-24],  This  algorithm  uses  predictor-corrector  updates  to  track  the  root  as  a 
function  of  z  and  r,  and  is  given  by 

4+1,+  =  4+{^(4,zky)rH-F(4,zky)-^(Ak2jyxzk+1  -zk) 

r)F 

-—(Ak2,zk,rk)(rk+1-rk)}, 

Ak+1  =max{0,A2+h+}. 

To  implement  the  algorithm,  we  take  advantage  of  the  known  functional  form  for  F  and  explic- 
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itly  compute  the  partial  derivatives, 


dX 
dX2 
8F 
8X2 
dX 
dr 
dF_ 
~8 r 
dX 
dz 
dF_ 
dz 


(P  +  X2s' T QS r 1  {- -s' T Q(S Xe(r)  - z)  -  S 1 QS x) , 

2(SX-z)tQ(S |^), 

OX2 

(P  +  X2SrQSyl  j-STeSQ}d2, 

T  dX 

2  (SX-z  +  r)TQ(S  —  +I3), 
or 

{P  +  X2STQSrlSTQSm2, 

T  dX 

2(5  X-z  +  r)T  Q(S  — —  /3), 


(12) 


where,  Xe(r)  =  Fir, 


and  It,  denotes  the  3x3  identity  matrix.  Note  that  S  FI  =  I2. 

Figure  3  illustrates  growth  distance  tracking.  For  the  first  20  iterations,  rk  is  held  constant  to 
enable  initial  convergence  of  the  algorithm.  Then,  rk  varies  through  the  virtual  net.  One  iteration 
of  the  Newton-Raphson  algorithm  per  value  of  rk  is  used  to  update  the  root,  X^+l.  Figure  3b 
demonstrates  that  the  growth  distance  tracking  is  accurate.  The  growth  distance  is  occasionally 
zero  indicating  an  overlap  between  several  r*  and  the  debris.  Figure  3c  illustrates  the  trajectory  of 
r*  in  three  dimensions. 

3.7  Thrust  Limit  on  Growth  Distance  Computations 

Suppose  that  the  thrust  limits  are  expressed  in  the  form  \\LU\\  <  1  for  an  appropriately  defined 
matrix  L  and  norm  ||  •  ||.  The  computational  procedures  to  determine  pu(r,K )  involve  solving  a 
bilevel  optimization  problem  where  \\L(KX  +  H(K)r)\\  is  maximized  subject  to  the  constraint  X  e 
C(r,K,a),  and  bisections  are  performed  on  the  value  of  a  so  that  the  maximum  value  is  driven  to 
1 .  As  we  demonstrate  in  this  section,  in  special  cases  this  computation  can  be  greatly  simplified. 

Suppose  that  the  thrust  constraints  are  prescribed  in  terms  of  polyhedral  norm  bounds,  specifi¬ 
cally 

r I  ( KX  +  Hr)  <  i/maxi  i  —  2,  •  •  • ,  m,  (13) 

where  c,-  are  the  vertices  of  the  unit  norm  polytope,  and  «max  is  the  norm  bound.  The  infinity  norm, 
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(C) 


Figure  3:  (a)  Components  of  r,  varying  versus  iteration  number,  (b)  Growth  distance  versus  itera¬ 
tion  number  computed  by  Newton-Raphson  algorithm,  (c)  The  trajectory  of  r  and  the  debris. 


for  instance,  has  m  =  6,  and 


1 

-1 

0 

e\  = 

0 

e2  = 

0 

£3  = 

1 

0 

0 

0 

0 

'  0  ' 

0 

*4  = 

-1 

^5  = 

0 

e6  = 

0 

0 

1 

-1 

In  the  case  of  non-polyhedral  norm  bounds,  such  as  the  2-norm,  an  approximation  by  a  polyhedral 
norm  bound  may  be  employed. 

The  thrust  limit  on  the  growth  distance  is  then  determined  based  on  solving,  for  i  =  1  ,•••,«,  the 
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optimization  problems 


maximize  ej (KX  +  Hr) 

subject  to  \(X-Xe(r))TP(X-Xe(r))<c.  (  ' 

If  the  value  of  c  is  found  for  which  the  solutions  X*  of  (15)  satisfy  ma Xj{eJ(KX*  +  Hr)}  =  nmax, 
then  pu(r,  K)  =  yfc. 

The  problem  (15)  can  be  solved  by  diagonalizing  P,  using  an  orthogonal  matrix,  V, 

P  =  VtAV,  A  =  diag  [A2,  ■■■,A2],Ai>  0. 

By  defining,  z  =  X-Xe(r),  and  £  so  that 


z= vTA-2<r, 


it  follows  that 


zTpz  =  £tA-12VPVtA~12£ 


ft- 


The  problem  (15)  can  now  be  re-written  as 


where 


maximize  hJ^  +  ejTr 
subject  to  C  <  c. 


hj  =  ejKVT  A- 2. 


(16) 


The  solution  to  the  constrained  maximization  problem  (16)  of  maximizing  the  inner-product  of 
two  vectors  over  a  unit  2-norm  ball  is  given  by 


&  = 

\M\ 


(17) 


where  ||  •  ||  denotes  the  vector  2-norm.  The  maximum  value  of  the  objective  function  in  (15)  is 
given  by 

\\hi\Wlc  +  ejTr. 


Consequently,  to  satisfy  (13),  we  let 


'  0, 

II 

.  1 
mm  - 

i  2 

(u  - ^TPr) 

Mmax  1  1 

\M  J 

if  3  i  :  Umax  <  ejYr , 


otherwise. 


(18) 
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Thus,  the  problem  of  finding  the  thrust  limit  on  the  growth  distance  for  polyhedral  norm  bounds 
has  an  explicit  solution  given  by  (18).  Even  though  the  computation  of  thrust  limits  on  the  growth 
distance  can  be  performed  offline  for  the  nominal  operating  conditions,  fast  computational  proce¬ 
dures  are  beneficial  in  case  of  thruster  failures,  degradations,  and  restrictions  on  thrust  directions 
(e.g.,  caused  by  the  presence  of  other  spacecraft  nearby),  all  of  which  can  lead  to  changing  con¬ 
straints  on  thrust  during  spacecraft  missions. 

We  note  that  the  condition  wmax  >  max/fc^Tr}  is  satisfied  if  the  available  thrust  can  maintain  the 
equilibrium  Xe(r )  in  steady-state.  We  also  note,  that,  based  on  the  form  of  T,  c  is  independent  of 
ry,  the  in-track  component  of  the  equilibrium  in  the  virtual  net.  Hence  the  computations  of  pu(r,K) 
need  only  be  performed  with  ry  =  0. 

When  a  spacecraft  does  not  have  independent  thrusters  in  x,  y  and  z  directions,  a  2-norm  thrust 
limit  is  more  practical.  Unfortunately,  (15)  is,  in  general,  a  non-convex  problem.  In  this  case,  the 
2-norm  bound  can  be  approximated  by  a  polyhedral  norm  bound  (13),  with  the  vertices  e,  selected 
on  the  unit  2-norm  ball  in  R 3.  We  note  that  higher  accuracy  of  this  approximation  requires  a  higher 
number  of  vertices  in  (13),  which  thus,  complicates  (18). 

Finally,  we  note  that  when  Av’s  are  treated  as  control  inputs,  the  thrust  limit  on  growth  distance 
is  induced  by  the  available  Av.  In  this  case,  computing  the  thrust  limited  growth  distance  is  com¬ 
pletely  analogous  to  computing  it  in  the  case  when  the  control  input  is  the  thrust  force  or  thrust 
acceleration. 


3.8  Connectivity  Graph  and  Graph  Search 

We  now  introduce  a  notion  of  connectivity  between  two  vertices  of  the  virtual  net,  r,-  e  J\f  and 
rj  e  J\f.  The  vertex  r,-  is  connected  to  the  vertex  rj  if  there  exists  a  gain  K  e  1C  such  that 


Xein)  e  intCirj, K,p*(rp  K,z)\  (19) 

where  int  denotes  the  interior  of  a  set.  The  connectivity  implies  that  a  spacecraft  located  close  to 
an  equilibrium  corresponding  to  r,-  can  transition  to  an  equilibrium  Xe(rj)  by  using  limited  thrust 
and  avoiding  collision  with  the  debris.  We  note  that  if  r,  is  connected  to  rj  this  does  not  imply 
that,  in  turn,  rj  is  connected  to  r(.  We  also  note  that  connectivity  depends  on  the  existence  of  an 
appropriate  control  gain  from  the  set  of  gains  1C  but  the  condition  (19)  does  not  need  to  hold  for  all 
gains. 

The  on-line  motion  planning  with  debris  avoidance  is  performed  according  to  the  following 
procedure  (for  simplicity,  described  here  for  the  case  of  a  single  debris): 

Step  1:  Determine  the  debris  location  and  shape  (i.e.,  z  and  Q ). 

Step  2:  By  using  fast  growth  distance  computations,  determine  the  thrust  limited  growth  dis¬ 

tance  based  on  (8),  with  pg  computed  online  and  pu  precomputed  offline. 
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Step  3:  Construct  a  graph  connectivity  matrix  between  all  r,-,  rj  e  Af.  In  the  graph  connectivity 

matrix,  if  two  vertices  are  not  connected,  the  corresponding  matrix  element  is  zero;  if  they 
are  connected  the  corresponding  matrix  element  is  1.  In  parallel,  build  the  control  gain 
selectivity  matrix,  which  identifies  the  index  of  the  highest  preference  gain  K  for  which  r,- 
and  rj  are  connected.  This  gain  will  be  applied  if  the  edge  connecting  r,-  and  rf  is  traversed. 

Step  4:  Perform  graph  search  to  determine  a  sequence  of  connected  vertices  r[k ]  e  Af  and  con¬ 

trol  gains  K[k]  e  1C,  k  =  1,  •  •  •  ,lp,  such  that  r[l]  satisfies  the  initial  constraints,  r[lp]  satisfies 
the  final  constraints,  and  the  path  length  lp  is  minimized. 

Per  the  above  algorithm,  graph  search  is  utilized  to  determine  the  minimum  number  of  equi¬ 
librium  hops  around  a  piece  of  debris.  After  the  path  has  been  determined  as  a  sequence  of  the 
set-points  and  the  corresponding  control  gains,  the  execution  of  the  path  proceeds  by  checking  if 
the  current  state,  X{t)  is  in  the  safe  positively  invariant  set  corresponding  to  the  next  reference  r+ 
and  next  control  gain  K+  in  the  sequence;  if  it  is,  then  the  controller  switches  to  this  reference  and 
control  gain: 


X(t)  G  C(r\K\p\r+,K+,z))  -»  r(t)  =  r+,  K(t )  =  K+ .  (20) 

3.9  Cost  Matrices 

As  described  in  the  previous  section,  the  connectivity  graph  matrix  is  comprised  of  ones  and  zeros, 
and  thus,  graph  search  results  in  a  minimum  length  path  between  desired  r,-,  rj  e  J\f . 

In  order  to  produce  time  efficient  and  thrust  efficient  paths,  offline  we  simulate  transitions 
between  all  r^rj  e  A f  for  each  K  e  1C  and  record  the  time  and  fuel  consumption  to  reach  a  box 
of  lm  around  the  target  vertex.  The  results  are  merged  into  time  and  fuel  matrices  that  store  the 
respective  minimum  value,  while  in  parallel,  the  control  selectivity  matrix  identifies  which  gain  K 
produced  said  minimum. 

Step  3  in  the  motion  planning  procedure  is  augmented  so  that  the  graph  connectivity  matrix 
is  multiplied  elementwise  with  a  desired  cost  matrix.  Vertices  that  are  not  connected  retain  a 
corresponding  matrix  element  of  zero,  while  vertices  that  are  connected  now  contain  a  matrix 
element  of  time  or  fuel  cost. 

3.10  Moving  Debris  Avoidance  Approach 

To  avoid  a  non- stationary  debris,  its  path  can  be  covered  by  a  union  of  a  finite  number  of  ellipsoidal 
sets, 

l=nd 

V={j0(zi,Qi),  (21) 

i=\ 

where  the  center  of  the  /th  set  is  denoted  by  zi  e  R3,  and  the  Ith  set  shape  is  defined  by  Qi  =  Qj  >  0. 
Then,  the  debris  avoidance  condition  for  the  closed-loop  trajectory  that  emanates  from  X(0)  with 
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the  set-point  r,-  and  gain  Kj  is  given  by 

X(0)  €  C(ru  Kj,p )  and  0(Zl,  Qi)  n  C(rf,  tf,,p)  =  0,  for  all  /  =  1 ,  •  •  • ,  nd.  (22) 

The  same  approach,  with  larger  nj,  can  be  used  to  handle  multiple  non-stationary  debris.  Note, 
however,  that  this  approach  is  conservative  as  it  does  not  account  for  the  debris  progressions  along 
their  paths  versus  time. 

Hence,  we  introduce  the  notion  of  time  into  the  problem;  whereas  a  transition  between  r,  and 
rj  might  not  be  feasible  at  time  t\ ,  based  on  the  motion  of  a  debris,  it  might  become  feasible  at 
time  t2 .  To  accommodate  moving  debris,  we  introduce  sets  Ck(r,K,p),  0  <  k  <  N,  defined  by  the 
following  relation, 


A(K)k(ck(r,  K,p )  -  {Xe(r)}j  c  ( C(r ,  K,p)  -  {Xe(r)}j,  (23) 

Note  that  if  X(0)  e  Ck(r,  K,p ),  then  X(l)  e  C*_i(r,  K,p),  X(2)  e  C*_2(r,  K,p ),  •  •  •  ,X(k)  e  C0(r,  K,p)  = 
C(r,K,p).  The  set  Ck(r,K,p )  can  be  much  larger  than  C(r,K,p );  any  states  in  Ck(r,K,p)  contract  to 
C(r,  K,p )  in  k  steps. 

We  now  define  connectivity  between  two  vertices  of  the  virtual  net,  r(-  e  J\f  and  rj  €  TV  at  a 
specified  time  to-  This  notion  is  based  on  the  fact  that  the  time  to  transition  from  any  state  in 
Cj^ir,  K,p)  to  C(r,  K,p)  is  less  or  equal  than  N  steps.  Suppose  that  the  debris  path  D(to  :  to  +  N  ■  H) 
has  been  predicted  over  the  N  ■  H  discrete  steps  from  the  time  instant  to,  where 

t=tr 

D(tk:tr)  =  [j0(z(t),Q(t)). 

t—tk 

The  node  r,  €  M  is  connected  to  r;  €  M  at  the  time  instant  tk  =  to  +  kN  if  there  exists  K  e  1C  such 
that 

D(tk  :  tk  +  N)C\  C(ri,  K,p)  =  0.  (24) 

The  node  r;  e  J\T  is  connected  to  node  rj  e  J\f  at  time  tk  if  there  exists  K  eK,  such  that 

D(tk  :tk  +  N)n  CN(rh  K,p )  =  0  (25) 

and 

C(ri,K,p)  c  CN(rj,K,p).  (26) 

The  connectivity  implies  that  a  spacecraft  located  close  to  an  equilibrium  corresponding  to  r,-, 
Xe(ri),  can  transition  close  to  an  equilibrium  Xe{rj)  between  the  time  instants  tk  and  tk  +  N  while 
avoiding  collision  with  the  debris.  We  note  that  if  r;  is  connected  to  rj  this  does  not  imply  that,  in 
turn,  rj  is  connected  to  r,-.  We  also  note  that  connectivity  depends  on  the  existence  of  an  appro¬ 
priate  control  gain  from  the  set  of  gains  K,  but  does  not  need  to  hold  for  all  gains.  Furthermore, 
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since  connectivity  depends  on  the  predicted  motion  of  the  debris,  connectivity/non-connectivity 
can  depend  on  time. 

The  on-line  motion  planning  with  debris  avoidance  is  performed  according  to  the  following 
procedure: 

Step  1:  Determine  the  debris  location,  shape  and  predict  the  debris  path  D(to  :  to  +  N  H) 

Step  2:  Construct  graph  connectivity  matrices  corresponding  to  tk,k  =  0, 1,  •  •  •  ,H.  In  the  graph 

connectivity  matrix,  if  two  vertices,  r,-  and  rj,  are  not  connected  at  tk,  the  corresponding 
matrix  element  is  zero;  if  they  are  connected  the  corresponding  matrix  element  is  1.  In 
parallel,  build  the  control  gain  selectivity  matrix,  which  identifies  the  index  of  the  highest 
preference  gain  K  for  which  4  and  rj  are  connected.  This  gain  will  be  applied  if  the  edge 
connecting  4  and  rj  is  traversed. 

Step  3:  Perform  graph  search  to  determine  a  sequence  r[tk]  €  Af  and  control  gains  K[k]  e  1C, 

k=  1,  •  •  • , lp,  such  that  r[t\]  satisfies  the  initial  constraints,  r[lp]  satisfies  the  final  constraints, 
and  the  path  length  lp  (or  another  cost  function  such  as  the  expected  fuel  consumption  or 
expected  maneuver  time)  is  minimized. 

Per  the  above  algorithm,  a  graph  search  is  utilized  to  determine  the  minimum  number  of  equi¬ 
librium  hops  around  a  debris  starting  at  to. 

Remark  1:  The  condition  (25)  is  conservative.  It  can  be  replaced  by  a  less  conservative 
condition, 

D(tk  :  4  +  m)  n  CN-m(rj,  K,p )  =  {0}, 
m  =  0, 1 ,  •  •  •  ,N, 

at  a  price  of  more  demanding  computations. 

Remark  2:  The  condition  (25)  is  checked  computationally  using  the  fast  growth  distance 
algorithm  described  in  Section  3.7.  The  intersection  is  empty  if  Cn  can  be  grown  before  it  touches 
D(tk  :  4  +  N).  This  fast  growth  distance  algorithm  is  essential  to  be  able  to  rapidly  construct  the 
connectivity  matrices. 

Remark  3:  In  our  simulations,  the  path  search  is  performed  using  the  standard  Dijkstra’s 
algorithm.  It  is  applied  to  a  lifted  graph  with  vertices  being  the  pairs  (4, 4). 

3.11  Bounded  Disturbances 

We  now  discuss  how  the  debris  avoidance  approach  can  be  extended  to  handle  bounded  distur¬ 
bances.  For  simplicity,  we  consider  the  case  of  multiple  stationary  debris.  Consider  the  system 

Xk+l  =AXk  +  BUk  +  Bw,  (27) 

where  w  eW,  W  is  the  convex  hull  of  wl  for  i  =  1, . . .  ,nw,  wl  are  the  vertices  of  a  disturbance  set, 
and  nw  is  the  number  of  vertices.  Note  that  W  is  a  compact  set. 
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The  positive  invariance  of  C(r,K,y),  y  >  0,  for  W  =  {0}  has  already  been  established.  When 
W  ±  {0},  it  can  be  shown  that  there  exists  ym in  such  that  the  set  is  positively  invariant  for  y  >  ym in. 
Note  that  ymin  =  ymin  (AT). 

Since  C(r, K,ymin(K))  is  disturbance  invariant,  it  contains  the  minimum  invariant  set  that  is  an 
attractor  for  closed-loop  trajectories,  as  long  as  r  and  K  are  maintained  at  constant  values.  Hence, 
in  the  case  of  bounded  disturbances,  connectivity  can  be  redefined  by  replacing  Xe(r,)  in  (19)  with 
C(r,,  A', ymm(K)).  Specifically,  the  vertex  r,-  e  N  is  connected  to  the  vertex  rj  e  J\f  if  there  exists 
K  e  1C  such  that 


C(ri,K,ym[n(K))  c  int C(r j, K,p* (r j, K,z)),  for  all  /  =  l,...,m.  (28) 

The  condition  (28)  ensures  that  a  switch  from  r,  to  rj  may  occur  and  that  subsequent  dynamics 
will  not  lead  to  collision  with  the  debris  once  X{t)  e  C(ri,K,ym{n(K)). 

To  compute  ym jn  under  all  possible  w  eW,  it  is  sufficient  to  examine  the  flow  at  the  vertices 
wl  of  the  disturbance  set  and  demonstrate  that  if  e  C(r,  K,  y{K))  and  w  e  {w\i  =  1 , . . . ,  nw},  then 
Xk+i  e  C(r,  K,y(K)).  The  value  ymjn  is  the  minimum  y  for  which  this  condition  holds. 

To  find  ymin  we  use  a  bilevel  optimization  strategy  where  the  inner  loop  solves  nw  optimization 
problems  numerically  with  respect  to  X, 

maximize  Ft(X)  =\(AX  +  BU  +  Bw')jF^-(AX  +  BU  +  Bw'), 

y]  (29) 

subject  to  \(X -  Xe(r))T P(K)(X  -  Xe(r))  <  yf, 

and  the  outer  loop  performs  bisections  on  each  y„  so  that  all  F/(X:Ty/)),  where  X*(ji)  denotes 
the  inner-loop  solution,  are  driven  to  1.  Thus,  ymin  =  min(y;)  for  i  =  1  ,...,nw.  Note  that  ymjn  is 
independent  of  equilibrium  r,  and  so  this  calculation  may  be  done  once  offline  for  each  K  eK.  and 
stored  onboard  for  real  time  implementation. 


4  RESULTS  AND  DISCUSSION 

Simulations  are  now  provided  that  illustrate  the  debris  avoidance  approach.  We  consider  a  nomi¬ 
nal  circular  orbit  of  850  km  and  discretize  the  HCW  equations  with  a  sampling  period,  AT,  of  30 
seconds.  We  construct  an  approximately  2  km  cubed  virtual  net.  We  let  K.  =  {Aj ,  A3,  where 
K\,K2,K?)  are  the  LQ  gains  associated  with  weight  matrices  Q  =  diag(100, 100, 100, 107, 107, 107), 
and  R\  =  2  x  105/3,  R2  =  2  x  107/3,  and  A3  =  2  x  109/3.  These  gains  are  chosen  to  represent  prefer¬ 
ences  for  fuel  considerations,  maneuver  time  considerations,  and  a  compromise  between  them.  We 
impose  a  maximum  thrust  constraint  of  10  N  in  each  axis.  In  all  simulations,  Dijkstra’s  algorithm 
is  used  to  find  the  shortest  cost  path  from  initial  node  to  final  node. 
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4.1  Static  Debris 


We  consider  an  ellipsoidal  set  0(zi,  Q\)  over-bounding  a  debris  centered  at  z\  =  [0.3  0.4  0.5]T  km, 
where  Q\  =  IOO/3.  We  use  the  technique  of  [25]  where  bisections  are  applied  to  (11)  to  determine 
the  growth  distance  to  the  debris  from  each  node  in  the  net.  The  spacecraft’s  initial  condition  is 
X(0)  =  W( /'ok  where  ro  =  [0.32  0  1.61]T  km.  The  target  equilibrium  node  is  X£,(0). 

Figure  4  shows  the  path  the  spacecraft  takes  under  closed-loop  control  in  order  to  avoid  the 
debris.  The  green  x  marks  the  initial  node.  The  blue  x  marks  the  final  node.  The  red  ellipsoid 
represents  the  debris.  The  blue  line  is  the  path  the  spacecraft  takes  in  order  to  avoid  the  debris.  The 
blue  ellipsoids  represent  the  invariant  sets  along  the  path.  The  spacecraft  is  able  to  complete  the 
desired  maneuver  well  within  maximum  thrust  constraints  while  successfully  avoiding  the  debris. 
In  Figure  5  we  rerun  the  simulation  for  a  grid  of  initial  conditions.  The  figure  clearly  demonstrates 
the  initial  conditions  for  which  the  maneuver  path  is  perturbed  from  that  which  the  spacecraft 
would  have  taken  had  there  been  no  debris. 


Figure  4:  (a)  Debris  avoidance  path  for  a  single  debris,  (b)  The  time  history  of  thrust  magnitude. 


Next,  we  add  a  second  debris  0(z2 ,  Qi)  centered  at  zi  =  [0.3  -  0.4  0.5]T,  where  Q2  =  IOO/3.  In 
calculating  the  growth  distance,  we  take  the  minimum  distance  to  each  of  0(z,i,  Qi),  i-  1,2.  Figure 
6  shows  the  path  the  spacecraft  takes  under  closed-loop  control  in  order  to  avoid  both  debris. 

4.2  Moving  Debris 

We  consider  the  case  of  a  non-stationary  debris  where  we  treat  its  motion  as  the  union  of  static 
debris  along  the  path  (21).  A  union  of  ellipsoidal  sets  over-bounds  the  debris’  motion,  where 
the  debris  positions  zi  are  generated  by  sampling  the  relative  motion  of  the  debris  with  the  initial 
condition  [0  0.5  0  0  0.0006  0]T,  and  where  Q,  =  2OO/3,  i=l ...  n(j.  The  spacecraft’s  initial  condition 
is  X(0)  =  Xe(ro),  where  ro  =  [0  1  0]T  km.  The  target  equilibrium  node  is  Xe(r(j),  where  ro  =  [0  - 
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Figure  5:  Debris  avoidance  paths  for  many  initial  conditions.  Each  green  x  marks  an  intial  condi¬ 
tion.  We  do  not  show  the  invariant  set  ellipsoids  for  visual  clarity. 


Figure  6:  (a)  Debris  avoidance  path  for  2  pieces  of  debris,  (b)  The  time  history  of  thrust  magnitude. 
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1  0]T  km.  We  use  the  single  gain  K2  and  do  not  include  fuel  or  time  cost  matrices  in  the  simulation, 
searching  for  a  minimum  length  path.  Figure  7  demonstrates  that  the  spacecraft  is  able  to  avoid 
the  closed  debris  path  by  ‘hopping’  under  it.  The  green  x  marks  the  initial  node.  The  blue  x  marks 
the  final  node.  The  red  ellipsoids  represent  the  debris  path.  The  blue  line  is  the  path  the  spacecraft 
takes  in  order  to  avoid  the  debris.  The  blue  ellipsoids  represent  the  maximally  grown  invariant  sets, 
C,  along  the  path. 

In  Figure  8,  we  repeat  the  simulation  for  time  efficient  and  thrust  efficient  paths  and  allow 
all  K  €  /C.  Table  1  summarizes  the  total  time,  thrust  and  nodes  traversed  for  the  three  paths. 
Note  that  the  minimum  length  path  now  ‘hops’  over  the  debris  path  instead  of  under  it,  as  now 
that  it  has  access  to  K\  it  finds  a  shorter  path.  Also  note  that  the  time  efficient  path  takes  longer  to 
complete  than  the  minimum  length  path.  While  the  cost  matrices  described  in  Section  3.9  calculate 
time  and  thrust  to  travel  between  all  vertices  in  the  virtual  net,  the  execution  of  the  path  does  not 
require  the  spacecraft  to  reach  intermediate  vertices,  rather,  switching  to  the  next  reference  once 
the  current  state  enters  the  next  reference’s  invariant  set  (20).  As  such,  the  cost  matrices  only 
provide  a  heuristic  for  selecting  efficient  paths.  In  Figure  9  we  require  the  paths  to  travel  through 
intermediate  vertices  to  show  that,  in  this  case,  the  cost  matrices  accurately  determine  efficient 
paths.  The  results  are  summarized  in  Table  2. 

Table  1:  Total  Time,  Thrust,  and  Nodes  Traversed  for  all  Maneuver  Paths  for  a  Union  of  Static 
Debris. 


Total  Time 

Total  Thrust 

Total  #  of  Nodes 

Gains  used 

Minimum  Length  Path 

2611.5  s 

1472.85  N-s 

6 

Ki 

Time  Efficient  Path 

2841  s 

1264.95  N-s 

6 

KuK2 

Thrust  Efficient  Path 

9177  s 

671.297  N-s 

11 

k2,k3 

Table  2:  Total  Time,  Thrust,  and  Nodes  Traversed  for  all  Maneuver  Paths  that  Travel  Through 
Intermediate  Nodes  for  a  Union  of  Static  Debris. 


Total  Time 

Total  Thrust 

Total  #  of  Nodes 

Gains  used 

Minimum  Length  Path 

10457.5  s 

3006.13  N-s 

6 

Ki 

Time  Efficient  Path 

9862  s 

2017.11  N-s 

6 

KuK2 

Thrust  Efficient  Path 

32812.5  s 

1083.58  N-s 

11 

k2,k3 
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Thrust,  N 


Figure  7:  (a)  Debris  avoidance  path  for  a  non-stationary  debris  using  the  union  method,  (b)  The 
time  history  of  thrust  magnitude,  (c)  Cumulative  thrust  vs  time. 
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*  Initial  Point 

*  Final  Point 
Shortest  Path 
Time  Eff  Path 
Fuel  Eff  Path 


(b) 


(c) 


Figure  8:  (a)  Multiple  debris  avoidance  paths  for  a  non- stationary  debris  using  the  union  method, 
(b)  The  time  history  of  thrust  magnitude,  (c)  Cumulative  thrust  vs  time. 
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*  Initial  Point 

*  Final  Point 
Shortest  Path 
Time  Eff  Path 
Fuel  Eff  Path 


Figure  9:  (a)  Debris  avoidance  paths  that  travel  via  intermediate  nodes  for  a  non- stationary  debris 
using  the  union  method,  (b)  The  time  history  of  thrust  magnitude,  (c)  Cumulative  thrust  vs  time. 
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We  now  repeat  the  simulations  taking  into  account  the  debris’  motion  as  a  function  of  time.  We 
use  the  single  gain  K2  and  do  not  include  fuel  or  time  cost  matrices  in  the  simulation,  searching 
for  a  minimum  length  path.  Figure  10  shows  that  the  graph  search  algorithm  is  able  to  find  a  path 
which  passes  through  the  debris’  path  but  avoids  collision  due  to  the  debris’  location  elsewhere 
at  the  specific  time  instant  at  which  the  spacecraft  path  crosses  the  debris’  path.  In  Figure  11, 
we  repeat  the  simulation  for  a  thrust  efficient  path.  Table  3  summarizes  the  total  time,  thrust  and 
nodes  traversed  for  the  two  paths.  Note  that  the  thrust  efficient  path  uses  more  thrust  than  the 
minimum  length  path.  In  Figure  12  we  require  the  paths  to  travel  through  intermediate  vertices  to 
show  that,  in  this  case,  the  thrust  cost  matrix  accurately  determines  an  efficient  path.  The  results 
are  summarized  in  Table  4. 

Table  3:  Total  Time,  Thrust,  and  Nodes  Traversed  for  all  Maneuver  Paths  using  the  Contractive 
Set  Approach. 


Total  Time 

Total  Thrust 

Total  #  of  Nodes 

Minimum  Length  Path 

4635.5  s 

4635.5  N-s 

7 

Thrust  Efficient  Path 

4703.5  s 

781.407  N-s 

7 

Table  4:  Total  Time,  Thrust,  and  Nodes  Traversed  for  all  Maneuver  Paths  that  Travel  Through 
Intermediate  Nodes  using  the  Contractive  Set  Approach. 


Total  Time 

Total  Thrust 

Total  #  of  Nodes 

Minimum  Length  Path 

13388.5  s 

2060.14  N-s 

7 

Thrust  Efficient  Path 

12657.5  s 

957.116  N-s 

7 

Finally,  we  run  the  simulation  for  the  case  of  bounded  disturbances.  We  consider  W  =  {w  : 
IMloo  <  s}  for  which  nw  =  8,  that  is,  disturbances  that  fit  in  a  box  of  magnitude  s.  In  Figure  13  we 
consider  a  uniform  distribution  of  disturbances,  for  s  =  0.1  N  and  s  =  0.2  N.  The  orange  ellipsoids 
represent  the  disturbance  invariant  sets  C(r,K,ymin(K)),  along  the  path.  The  spacecraft  is  able  to 
safely  avoid  the  debris’  path  despite  being  subjected  to  disturbances. 

5  CONCLUSIONS 

We  described  a  technique  for  spacecraft  maneuver  planning  that  uses  positively-invariant  sets  in  or¬ 
der  to  avoid  collisions  with  debris,  while  adhering  to  specified  thrust  limits.  The  approach  is  based 
on  hopping  between  neighborhoods  of  equilibria  in  a  virtual  net,  and  maintaining  the  spacecraft 
trajectory  within  a  tube  formed  by  safe  positively-invariant  sets.  For  the  case  where  thrust  limits 
can  be  specified  as  polyhedral  norm  bounds,  we  have  shown  that  the  thrust  limit  on  the  growth 
distance  can  be  easily  computed;  it  is,  in  fact,  feasible  to  perform  these  computations  onboard  a 
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Figure  10:  (a)  Debris  avoidance  path  for  a  non- stationary  debris  using  the  contractive  set  approach, 
(b)  The  time  history  of  thrust  magnitude,  (c)  Cumulative  thrust  vs  time. 
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Figure  11:  (a)  Multiple  debris  avoidance  paths  for  a  non- stationary  debris  using  the  contractive  set 
approach,  (b)  The  time  history  of  thrust  magnitude,  (c)  Cumulative  thrust  vs  time. 
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Figure  12:  (a)  Debris  avoidance  paths  that  travel  via  intermediate  nodes  for  a  non- stationary  debris 
using  the  contractive  set  approach,  (b)  The  time  history  of  thrust  magnitude,  (c)  Cumulative  thrust 
vs  time. 
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Thrust,  N 


Figure  13:  (a)  Debris  avoidance  path  for  a  non- stationary  debris  under  uniform  random  disturbance 
with  s  =  0.1  N.  (b)  s  =  0.2  N.  (c),  (d)  Time  histories  of  thrust  magnitude. 
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spacecraft  in  order  to  account  for  thruster  failure  or  degradation.  We  described  an  extension  in 
the  presence  of  moving  debris  using  contractive  constraint  admissible  sets  in  order  to  avoid  col¬ 
lisions.  Lastly,  we  illustrated  that  the  approach  can  be  extended  to  include  unmeasured  bounded 
disturbances. 

Future  Work 

Developing  cost  matrices  that  accurately  capture  the  cost  of  maneuvers  that  do  not  travel  through 
intermediate  nodes  is  a  topic  for  future  investigation. 

The  constraint-admissible  positively  invariant  set  method  can  be  extended  to  attitude  control 
on  SO(3)  that  is  capable  of  handling  inequality  constraints  associated  with  control  authority  limits 
and  conical  keep-out  zones.  The  controller  would  use  a  supervisory  strategy  with  an  inner-loop 
Lyapunov  SO(3)-based  controller,  such  as  the  inertia-free  controllers  presented  in  this  work,  and 
an  outer  loop  set-point  guidance  based  on  positively  invariant  constraint  admissible  sets  with  real¬ 
time  graph  search.  The  combined  methodology  would  reduce  the  search  space  of  possible  attitude 
maneuver  solutions  and  effectively  handle  constraints. 

Future  work  will  also  consider  ways  to  apply  the  safe  positively  invariant  set  method  to  non¬ 
spacecraft  problems,  such  as  ground  and  other  autonomous  vehicles. 
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APPENDIX 


This  appendix  provides  the  necessary  background  on  spacecraft  relative  motion  orbital  dynamics. 

A.l  Spacecraft  Relative  Motion  Orbital  Dynamics 

In  traditional  relative  motion  problems,  an  approaching  spacecraft  is  maneuvered  close  to  a  target 
spacecraft  in  a  nominal  orbit.  The  target  spacecraft  is  assumed  to  be  at  the  origin  of  Hill’s  frame. 
See  Fig  A.  1 . 


Figure  A.  1 :  Hill’s  frame. 


A.1.1  Nonlinear  equations  of  motion 

The  relative  position  vector  of  the  spacecraft  with  respect  to  a  target  location  on  an  orbit  is  ex¬ 
pressed  as 

8r  =  xi+yj  +  zk, 

where  x,  y  and  z  are  the  components  of  the  position  vector  of  the  spacecraft  relative  to  the  target 
location  and  i,  j,  k  are  the  unit  vectors  of  the  Hill’s  frame.  The  Hill’s  frame  has  its  x-axis  along 
the  orbital  radius,  y-axis  orthogonal  to  the  x-axis  and  in  the  orbital  plane,  and  z-axis  orthogonal  to 
orbital  plane. 

The  position  vector  of  the  spacecraft  with  respect  to  the  center  of  the  Earth  is  given  by  R  = 

Ro  +  Sr,  where  Rq  is  the  nominal  orbital  position  vector.  The  nonlinear  equation  of  motion  for  the 
spacecraft  (relative  to  an  inertial  frame)  is  given  by 

R  =  -^  +  —F,  (A.l) 

R i  mc 
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where  F  is  the  vector  of  external  forces  applied  to  the  spacecraft,  R  =  \R\,  mc  is  the  mass  of  the 
spacecraft,  and  pi  is  the  gravitational  constant. 


A.1.2  Linearized  equations  on  circular  orbits 


For  5r  «  R,  the  linearized  CWH  equations  [2]  approximate  the  relative  motion  of  the  spacecraft 
on  a  circular  orbit  as 

9  Fx 

x-3nzx-2ny  =  — , 
mc 


y  +  2  nx 


in 


(A.2) 


z  +  n  z  = 


mr 


where  Fx,Fy,Fz  are  components  of  the  external  force  vector  (excluding  gravity)  acting  on  the 

spacecraft,  and  n  =  /5  denotes  the  mean  motion  of  the  nominal  orbit.  The  linearized  dynamics 

\  Ro 

account  for  differences  in  gravity  between  the  spacecraft  and  nominal  orbital  location,  and  for 
relative  motion  effects.  The  spacecraft  relative  motion  dynamics  in  the  orbital  plane  (x  and  y) 
and  in  the  out-of-orbital  plane  (z)  are  decoupled.  The  in-plane  dynamics  are  Lyapunov  unstable 
(2  eigenvalues  at  the  origin  and  2  eigenvalues  on  the  imaginary  axis  at  n ),  while  the  out-of-plane 
dynamics  are  Lyapunov  stable  (2  eigenvalues  on  the  imaginary  axis  at  n ).  The  in-plane  dynamics 
are  completely  controllable  from  Fy  input  but  are  not  controllable  from  Fx  input.  The  out-of-plane 
dynamics  are  controllable  from  Fz  input. 

Assuming  a  sampling  period  of  AT  sec,  we  can  convert  the  model  (A.2)  to  a  discrete-time  form 


Xk+l=AXk  +  BUk,  (A3) 

where  Xk  =  [xk,  yk,  Zk,  At,  jk ,Zk]T  is  the  state  at  time  step  k  e  Z+,  Uk  =  [Fxy,  Fyuk,  FZ^]T  is  the 

rAT 

control  vector  of  thrust  forces  at  the  time  step  k  e  Z+,  and  A  =  exp(AcAT),  B  =  J()  exp(Ac(AT  - 
r))drBc  are  the  discretized  matrices  obtained  based  on  the  continuous-time  system  realization 
(Ac,  Bc)  in  (A.2).  Alternatively,  the  control  vector  U  can  represent  an  instantaneous  change  in  the 
velocity  of  the  spacecraft,  Av,  induced  by  thrust,  with  an  appropriately  re-defined  5-matrix, 

0  0 
0  0 
0  0 
0  0  ' 

1  0 
0  1 


5av  =  e 


AnAT 


0 

0 

0 

1 

0 
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A.1.3  Linearized  equations  on  elliptic  orbits 

For  generic  elliptic  orbits  of  arbitrary  eccentricity,  the  linearization  of  these  equations  is  described 
by  linear  time-varying  equations  [26], 


Fx 

—  =  6x- 
mc 


2 11  h2  )s  2(vo(t)-  Ro(t))h 

+  — —  \Sx  + - 7 - Sy 


yRp)  Rp ) 
h 

2  ~^—6y, 

Rp) 


Rp) 


Fy 

—  =  Sy  + 
mc 


+  2- 


F 


h 2 


yRp)  Rp), 
h 

-ox. 


c  2(vo(t)-Ro(i))h 
oy - - - ox 

Rp) 


(A.4) 


Rp) 


tz  x"  a.  F 

—  =6z+  —7 — 
mc  Rp) 


Sz, 


where  5x,  8y  and  Sz  are  (relative)  coordinates  of  the  spacecraft  in  Hill’s  frame,  Fx,Fy,Fz  are  com¬ 
ponents  of  the  external  force  vector  (excluding  gravity)  acting  on  the  spacecraft,  h  is  the  orbit 
angular  momentum,  Ro(t)  is  the  nominal  time-varying  orbital  radius,  and  vo(0  is  the  nominal  time- 
varying  orbital  velocity.  Equation  (A. 2)  assumes  that  the  target  spacecraft  motion  is  in  an  ideal 
Keplerian  orbit;  if  its  motion  is  affected  by  perturbations,  Fx,  Fy,  Fz  have  to  be  modified  to  account 
for  these  perturbations  [27].  We  assume  that  Fx,Fy,Fz  are  thrust  forces  that  can  be  realized  via 
on-board  thruster  on-off  time  allocation  and  attitude  control  system  commands  [27]. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


CWH  Clohessy-Wiltshire-Hill 
LTV  linear  time- varying 
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